input data
source(file = '../inputProcessing/date_input.r')
rep <- 1e2
# rep_sim <- 1e4
for (date_num in 1:length(dates_input)){
date_week_finishing <- as.Date(dates_input[date_num],format = '%Y-%m-%d')
output <- list()
# date_week_finishing <- as.Date('29/03/2020',format = '%d/%m/%Y')
# # si <- 2
# rep <- 1e2
# rep_sim <- 1e2
for (Mdata in c('A','G')){
inputs<- readRDS(file=paste0(
'../../Mobility.2.0.Save/Saved_file/inputProcessing/Rdata/',
date_week_finishing,'_inputs_from_',Mdata,'.rds'))
for (si in 1:2){
########################################
####inputs
D <- inputs$D
M <- inputs$M
Ot <- inputs$Ot[[si]]$Ot
W <- inputs$W[[si]]$W
H <- inputs$H
SI <- inputs$SI[[si]]
delta_id <- inputs$delta_id
#
country <- names(D)[-1]
mD <- as.matrix(D[,-1])
N_geo <- ncol(D)-1
N_days <- nrow(D)
N_week <- nrow(D)/7
# prerp MCMC
# epi
R0 <- 5
# for risk
# b <- 0
o <- 1
theta0 <- c(rep(R0,N_geo*N_week),
rep(o,N_geo*N_week))#,
# rep(b,N_geo))
# theta0[2] <- 2
# theta0[3] <- 3
prior_theta <- matrix(c(rep(c(0,5),N_geo*N_week),rep(c(0,1e5),N_geo*N_week)),
length(theta0),2, byrow=TRUE)
# parameter names
f0 <- function(x) paste0('R',x)
f1 <- function(x,y) paste0(x,'_',y)
# f1 <- function(x) paste0('beta_',x)
# f2 <- function(x) paste0('Over_',x)
n_t<- rep(sapply(1:N_week,f0),N_geo)#, sapply(country,f1))#, sapply(country,f2))
n_t<- rep(country,each = N_week)#, sapply(country,f1))#, sapply(country,f2))
n_t<-c(n_t,'over')
# sd dev for proposal
sigma <- rep(1e-1,length(theta0))
########################################
# useful functions
sapply(paste0('Rscript/',(list.files('Rscript/'))),FUN = source)
# R_mat <- matrix(rep(theta0,each=7),nrow = N_days, ncol = N_geo,byrow = FALSE)
########################################
#run MCMC
# res <- MCMC_iter(iter = rep, theta0 = theta0, s = sigma)
res <- MCMC_full(iter = rep,
theta0 = theta0,
s = sigma,
repli_adapt = 10,
within_iter = rep/10)
########################################
## check convergence
print('########################################')
print(paste0('check convergence',date_week_finishing,' - ',Mdata,' - ',si))
Acc <- colSums(diff(res$theta)!=0)/(rep-1)
# Acc
if(rep >1e2) {
f <- round(seq(1,rep,length.out = 1e2))
}else{
f <- 1:rep
}
plot(res$logL[,1],main=paste0('DIC ',res$DIC[1],', P ',res$DIC[2]))
layout(matrix(1:4,2,2))
for (i in 1:(length(theta0)/2+1)){
if(i == (length(theta0)/2+1)){
plot(res$theta[f,i],
main = paste0(n_t[i],' - ',round(Acc[i]*100)) )
}else{
plot(res$theta[f,i],
main = paste0(n_t[i],' - ',round(Acc[i]*100)) , ylim = prior_theta[i,])
}
}
# apply(res$theta,2,quantile,c(.5,.025,.975))
########################################
## Baseline model: estimate Rt
summary <- apply(res$theta,2,quantile,c(.5,.025,.975))
median <- matrix(rep(summary[1,],each=7),nrow = N_days, ncol = N_geo,byrow = FALSE)
low <- matrix(rep(summary[2,],each=7),nrow = N_days, ncol = N_geo,byrow = FALSE)
up <- matrix(rep(summary[3,],each=7),nrow = N_days, ncol = N_geo,byrow = FALSE)
mean_r <- matrix(rep(apply(res$theta,2,mean),each=7),nrow = N_days, ncol = N_geo,byrow = FALSE)
sd_r <- matrix(rep(apply(res$theta,2,sd),each=7),nrow = N_days, ncol = N_geo,byrow = FALSE)
# format outputs
temp <- D
temp[,-1] <- NA
results_baseline <- list(median_R = temp,
low_R = temp,
up_R = temp,
CV = temp,
CV_below.2 = temp,
mean = temp)
for (i in 1:N_geo){
f <- which(cumsum(D[,i+1])>0)
results_baseline$median_R[f,i+1] <- median[f,i]
results_baseline$low_R[f,i+1] <- low[f,i]
results_baseline$up_R[f,i+1] <- up[f,i]
results_baseline$CV[f,i+1] <- sd_r[f,i]/mean_r[f,i]
results_baseline$CV_below.2[f,i+1] <- results_baseline$CV[f,i+1] < .2
results_baseline$mean[f,i+1] <- mean_r[f,i]
# results_baseline$M_D[,i+1] <- M_D[,i]
}
results_baseline$DIC <- res$DIC
output[[si]] <- list(resEpi = results_baseline)
}
saveRDS(object = output,
file = paste0('../../Mobility.2.0.Save/Saved_file/EpiEstim/Rdata/',
date_week_finishing,'_output_from_',Mdata,'.rds' ))
}
}
## [1] "########################################"
## [1] "check convergence2020-05-03 - A - 1"









































































































































































## [1] "########################################"
## [1] "check convergence2020-05-03 - A - 2"









































































































































































## [1] "########################################"
## [1] "check convergence2020-05-03 - G - 1"































































































































## [1] "########################################"
## [1] "check convergence2020-05-03 - G - 2"































































































































## [1] "########################################"
## [1] "check convergence2020-05-10 - A - 1"











































































































































































## [1] "########################################"
## [1] "check convergence2020-05-10 - A - 2"











































































































































































## [1] "########################################"
## [1] "check convergence2020-05-10 - G - 1"























































































































































## [1] "########################################"
## [1] "check convergence2020-05-10 - G - 2"























































































































































